Disclaimer

This work is very preliminary as I get back into the coding swing of things. Data wrangling and figure generation will be done via R, but the rest of the project will be done using good ol’ microsoft products. This is just an entry point into data crunching and should by no means be considered a final product.

3/29 update:

slight error in methodology- should be corrected now. Git not showing up (but maybe it will after this next push!)

All COOP

Coop Data

NOAA/NWS Cooperative Observer Network

For the Yampa area:

STEAMBOAT SPRINGS CO7936 HAYDEN CO3867 WALDEN CO8756 KREMMLING CO4664 CRAIG 4 SW CO1932

For the Rio Grande area:

DEL NORTE 2E CO2184 HERMIT 7 ESE CO3951 PAGOSA SPRINGS CO6258 SILVERTON CO7656 VALLECITO DAM CO8582

Downloading from the Iowa Environmental Mesonet.

COOP_stations <- read.csv("C:/Users/13074/Documents/ESS580/thesis_project/All_COOP_stations/data_raw/COOP_stations.csv", header = TRUE)

#str(COOP_stations)
COOP_stations$Date <- ymd(COOP_stations$day)

COOP_stations_clean <- COOP_stations %>% # filter for the timeframe
  addWaterYear() %>%
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days")))) %>%   na.omit()

COOP_stations_clean <- COOP_stations_clean %>% 
  mutate(avg_T_c = (highc+lowc)/2) %>% 
  filter(waterYear <= 2022)


write.csv(COOP_stations_clean,"C:/Users/13074/Documents/ESS580/thesis_project/All_COOP_stations/data_clean/COOP_stations_clean.csv", row.names = FALSE)

ggplot(COOP_stations_clean, aes(x = Date, y = avg_T_c, color = station_name)) +
  geom_line() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

# Just checking.

Station methodology

Two stations (Del Norte and Steamboat) will be analyzed for longer spans of time (1900 - 2022), while the remaining stations will be analyzed over the same period of record that the SNOTEL stations operate at (~1980s -2022). Most SNOTEL stations analyzed do not have reliable temperature data before WY 1987 & 1984 in the northern and southern areas respectively. For this purpose we will look at data after WY 1984.

Standard Deviation

To figure out the standard deviation for each year, I want the “residual” for each daily value.

The standard deviation will be the daily residual minus the mean of the residuals by water year, summed and squared, then divided by the number of observations minus one. The square root of the resulting value of which is thus the standard deviation for the water year.

Northern Stations

HAYDEN CO3867

Detrending Data

#average water year temperature

CO3867_yearly_wy_aver <- COOP_stations_clean %>%
  filter(station == "CO3867") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO3867_daily_wy_aver <- CO3867_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

CO3867_daily_wy_aver <- CO3867_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO3867_daily_wy_aver$aver_day_temp))

#str(CO3867_daily_wy_aver)
# try to show all years as means. 
CO3867_daily_wy_aver2 <- CO3867_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO3867_daily_wy_aver2$date_temp <- signif(CO3867_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO3867_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO3867_standard_dev <- CO3867_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

#checking the detrend

detrend <- CO3867_standard_dev %>% 
  filter(waterYear == 1984)

ggplot(detrend, aes(x = waterDay, y = residual))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

mean(CO3867_standard_dev$residual)
## [1] -2.129571e-16

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO3867_standard_dev_all <- CO3867_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO3867_standard_dev_all <- CO3867_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3867_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.022028
1985 4.332474
1986 4.259993
1987 3.864474
1988 3.872949
1989 4.430677
1990 3.467251
1991 4.070296
1992 3.897746
1993 3.885791
1994 3.705640
1995 4.075178
1996 3.880133
1997 4.302880
1998 3.625407
1999 4.119806
2000 3.653395
2001 3.869937
2002 3.858667
2003 3.810870
2004 3.825051
2005 3.738393
2006 4.139872
2007 4.108651
2008 3.803621
2009 3.796570
2010 3.884209
2011 4.101922
2012 3.393176
2013 4.428168
2014 4.037080
2015 4.289606
2016 3.782591
2017 4.289811
2018 3.584915
2019 3.899139
2020 4.167912
2021 3.907259
2022 3.736448 
#CALLING THIS something different
CO3867_all_V2 <- ggplot(CO3867_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')


print(CO3867_all_V2)

CO3867_sd_mk <- mk.test(CO3867_standard_dev_all$sd_2)
print(CO3867_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO3867_standard_dev_all$sd_2
## z = -0.62904, n = 39, p-value = 0.5293
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -53.00000000 6833.66666667   -0.07152497
CO3867_sd_sens <- sens.slope(CO3867_standard_dev_all$sd_2)
print(CO3867_sd_sens)
## 
##  Sen's slope
## 
## data:  CO3867_standard_dev_all$sd_2
## z = -0.62904, n = 39, p-value = 0.5293
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.009171542  0.004831871
## sample estimates:
##  Sen's slope 
## -0.002628673

Summer temperature standard deviation

CO3867_standard_dev_all_summer <- CO3867_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO3867_standard_dev_all_summer <- CO3867_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3867_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.445506
1985 2.487115
1986 2.093323
1987 2.064644
1988 1.999852
1989 2.201513
1990 2.354305
1991 1.726839
1992 2.445651
1993 2.402066
1994 2.187613
1995 2.742999
1996 1.830182
1997 2.114356
1998 2.647702
1999 1.593625
2000 2.006863
2001 2.305292
2002 2.337941
2003 2.508134
2004 2.411519
2005 2.542864
2006 1.912627
2007 2.180150
2008 2.573272
2009 2.052403
2010 2.152240
2011 1.879349
2012 1.832650
2013 1.858680
2014 1.925874
2015 2.319542
2016 2.261139
2017 1.992751
2018 2.142957
2019 2.431490
2020 2.439636
2021 2.491377
2022 2.135380 
ggplot(CO3867_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3867 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO3867_sd_mk_summer <- mk.test(CO3867_standard_dev_all_summer$sd_2)
print(CO3867_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO3867_standard_dev_all_summer$sd_2
## z = -0.26613, n = 39, p-value = 0.7901
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -23.00000000 6833.66666667   -0.03103914
CO3867_sd_sens_summer <- sens.slope(CO3867_standard_dev_all_summer$sd_2)
print(CO3867_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO3867_standard_dev_all_summer$sd_2
## z = -0.26613, n = 39, p-value = 0.7901
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.009457585  0.007800624
## sample estimates:
##  Sen's slope 
## -0.001356549

Winter temperature standard deviation

CO3867_standard_dev_all_winter <- CO3867_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO3867_standard_dev_all_winter <- CO3867_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO3867_standard_dev_all_winter <- CO3867_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3867_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.635204
1985 5.182071
1986 5.469629
1987 4.796823
1988 4.836877
1989 5.823914
1990 4.130722
1991 5.328556
1992 4.540552
1993 4.932496
1994 4.658666
1995 5.138381
1996 4.795413
1997 5.170437
1998 4.381221
1999 5.008785
2000 4.202000
2001 4.604372
2002 4.778003
2003 4.551116
2004 4.517692
2005 4.535750
2006 5.336493
2007 5.303816
2008 4.398136
2009 4.636602
2010 4.714850
2011 5.481038
2012 4.151779
2013 5.469158
2014 5.025790
2015 5.586965
2016 4.556407
2017 5.268566
2018 4.381153
2019 4.739618
2020 4.139964
2021 4.082859
2022 4.503057 
ggplot(CO3867_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3867 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO3867_sd_mk_winter <- mk.test(CO3867_standard_dev_all_winter$sd_2)
print(CO3867_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO3867_standard_dev_all_winter$sd_2
## z = -1.5242, n = 39, p-value = 0.1275
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##          S       varS        tau 
## -127.00000 6833.66667   -0.17139
CO3867_sd_sens_winter <- sens.slope(CO3867_standard_dev_all_winter$sd_2)
print(CO3867_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO3867_standard_dev_all_winter$sd_2
## z = -1.5242, n = 39, p-value = 0.1275
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.024252074  0.003511472
## sample estimates:
##  Sen's slope 
## -0.009536879

Spring and Fall temperature standard deviation

CO3867_standard_dev_all_spring <- CO3867_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO3867_standard_dev_all_spring <- CO3867_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3867_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.356153
1985 3.470221
1986 3.491535
1987 3.718489
1988 3.623295
1989 3.648456
1990 3.275793
1991 3.705977
1992 2.900631
1993 3.143210
1994 3.380657
1995 2.795633
1996 3.661786
1997 4.237358
1998 2.682349
1999 3.918743
2000 3.616376
2001 3.400696
2002 3.042341
2003 4.204434
2004 3.376469
2005 3.381613
2006 3.221596
2007 3.496825
2008 3.626653
2009 3.331802
2010 3.794990
2011 2.986717
2012 3.502668
2013 4.731100
2014 4.004636
2015 3.111478
2016 3.413607
2017 3.976452
2018 2.941087
2019 3.536352
2020 4.120725
2021 3.690872
2022 3.738474 
ggplot(CO3867_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3867 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO3867_sd_mk_spring <- mk.test(CO3867_standard_dev_all_spring$sd_2)
print(CO3867_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO3867_standard_dev_all_spring$sd_2
## z = 0.53226, n = 39, p-value = 0.5945
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 4.500000e+01 6.833667e+03 6.072874e-02
CO3867_sd_sens_spring <- sens.slope(CO3867_standard_dev_all_spring$sd_2)
print(CO3867_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO3867_standard_dev_all_spring$sd_2
## z = 0.53226, n = 39, p-value = 0.5945
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.009394674  0.016109893
## sample estimates:
## Sen's slope 
## 0.002947839

Fall temperature standard deviation

CO3867_standard_dev_all_fall <- CO3867_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO3867_standard_dev_all_fall <- CO3867_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO3867_standard_dev_all_fall <- CO3867_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3867_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.032027
1985 4.150946
1986 3.434685
1987 2.691224
1988 3.030391
1989 2.733709
1990 3.052697
1991 3.249156
1992 3.863224
1993 2.966372
1994 2.763922
1995 3.015554
1996 3.607809
1997 4.439213
1998 3.429782
1999 3.144617
2000 3.501683
2001 3.013489
2002 3.080980
2003 2.322525
2004 3.707258
2005 2.212636
2006 3.561168
2007 3.023899
2008 3.535942
2009 3.059430
2010 3.662327
2011 2.632570
2012 2.800936
2013 3.502222
2014 3.703235
2015 3.030694
2016 3.022735
2017 3.813281
2018 3.318303
2019 3.439630
2020 5.467131
2021 4.905616
2022 3.150385 
ggplot(CO3867_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3867 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO3867_sd_mk_fall <- mk.test(CO3867_standard_dev_all_fall$sd_2)
print(CO3867_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO3867_standard_dev_all_fall$sd_2
## z = 1.379, n = 39, p-value = 0.1679
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  115.0000000 6833.6666667    0.1551957
CO3867_sd_sens_fall <- sens.slope(CO3867_standard_dev_all_fall$sd_2)
print(CO3867_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO3867_standard_dev_all_fall$sd_2
## z = 1.379, n = 39, p-value = 0.1679
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.003868177  0.025543235
## sample estimates:
## Sen's slope 
## 0.009597072

WALDEN CO8756

Detrending Data

#average water year temperature

CO8756_yearly_wy_aver <- COOP_stations_clean %>% 
  filter(station == "CO8756") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO8756_daily_wy_aver <- CO8756_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

CO8756_daily_wy_aver <- CO8756_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO8756_daily_wy_aver$aver_day_temp))

#str(CO8756_daily_wy_aver)
# try to show all years as means. 
CO8756_daily_wy_aver2 <- CO8756_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO8756_daily_wy_aver2$date_temp <- signif(CO8756_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO8756_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO8756_standard_dev <- CO8756_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO8756_standard_dev$residual)
## [1] 1.886111e-16

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO8756_standard_dev_all <- CO8756_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO8756_standard_dev_all <- CO8756_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8756_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.014145
1985 4.079450
1986 4.282738
1987 3.811980
1988 4.225936
1989 4.494618
1990 3.759698
1991 4.029922
1992 3.892915
1993 3.809733
1994 3.813734
1995 4.043696
1996 4.254210
1997 4.315567
1998 3.698211
1999 4.165986
2000 3.638490
2001 3.907310
2002 3.720634
2003 4.069985
2004 4.128364
2005 4.037234
2006 4.693032
2007 4.870152
2008 4.485365
2009 4.238196
2010 4.835066
2011 4.508581
2012 3.950096
2013 5.182045
2014 4.306671
2015 5.006851
2016 4.510417
2017 4.863560
2018 4.205117
2019 4.533785
2020 4.964315
2021 4.468842
2022 5.071546 
CO8756_all <- ggplot(CO8756_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8756_all

CO8756_sd_mk <- mk.test(CO8756_standard_dev_all$sd_2)
print(CO8756_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO8756_standard_dev_all$sd_2
## z = 3.7984, n = 39, p-value = 0.0001456
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  315.0000000 6833.6666667    0.4251012
CO8756_sd_sens <- sens.slope(CO8756_standard_dev_all$sd_2)
print(CO8756_sd_sens)
## 
##  Sen's slope
## 
## data:  CO8756_standard_dev_all$sd_2
## z = 3.7984, n = 39, p-value = 0.0001456
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.01215542 0.03248087
## sample estimates:
## Sen's slope 
##   0.0222248

Summer temperature standard deviation

CO8756_standard_dev_all_summer <- CO8756_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO8756_standard_dev_all_summer <- CO8756_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8756_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.071927
1985 2.146719
1986 1.806037
1987 2.285426
1988 2.021938
1989 1.971525
1990 2.294354
1991 1.587187
1992 2.264433
1993 2.245478
1994 2.061340
1995 2.591182
1996 1.699686
1997 1.948128
1998 2.661193
1999 1.635943
2000 2.027596
2001 2.377058
2002 2.099655
2003 2.296450
2004 2.068967
2005 2.514033
2006 2.271628
2007 2.348102
2008 2.463497
2009 2.140242
2010 2.352250
2011 1.883893
2012 2.139116
2013 2.101973
2014 2.077100
2015 2.399835
2016 2.286265
2017 2.171841
2018 2.427565
2019 2.358237
2020 2.339032
2021 2.541649
2022 2.534067 
ggplot(CO8756_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8756 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO8756_sd_mk_summer <- mk.test(CO8756_standard_dev_all_summer$sd_2)
print(CO8756_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO8756_standard_dev_all_summer$sd_2
## z = 2.7581, n = 39, p-value = 0.005814
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  229.0000000 6833.6666667    0.3090418
CO8756_sd_sens_summer <- sens.slope(CO8756_standard_dev_all_summer$sd_2)
print(CO8756_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO8756_standard_dev_all_summer$sd_2
## z = 2.7581, n = 39, p-value = 0.005814
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.002732592 0.016007603
## sample estimates:
## Sen's slope 
## 0.008485008

Winter temperature standard deviation

CO8756_standard_dev_all_winter <- CO8756_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO8756_standard_dev_all_winter <- CO8756_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO8756_standard_dev_all_winter <- CO8756_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8756_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.956217
1985 5.038064
1986 5.854781
1987 4.837625
1988 5.372973
1989 6.069938
1990 4.774398
1991 5.114737
1992 4.634544
1993 4.942967
1994 4.990211
1995 5.058208
1996 5.584260
1997 5.498957
1998 4.688672
1999 5.196081
2000 4.347414
2001 4.811708
2002 4.765572
2003 5.315767
2004 5.201313
2005 5.206506
2006 6.176262
2007 6.480296
2008 5.587113
2009 5.015334
2010 5.982714
2011 6.065248
2012 4.790433
2013 6.635796
2014 5.346858
2015 6.781144
2016 5.881606
2017 6.256835
2018 5.112722
2019 5.679938
2020 5.523989
2021 4.883252
2022 6.623482 
ggplot(CO8756_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8756 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO8756_sd_mk_winter <- mk.test(CO8756_standard_dev_all_winter$sd_2)
print(CO8756_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO8756_standard_dev_all_winter$sd_2
## z = 2.5887, n = 39, p-value = 0.009633
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  215.0000000 6833.6666667    0.2901484
CO8756_sd_sens_winter <- sens.slope(CO8756_standard_dev_all_winter$sd_2)
print(CO8756_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO8756_standard_dev_all_winter$sd_2
## z = 2.5887, n = 39, p-value = 0.009633
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.00619369 0.04107522
## sample estimates:
## Sen's slope 
##  0.02169687

Spring and Fall temperature standard deviation

CO8756_standard_dev_all_spring <- CO8756_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO8756_standard_dev_all_spring <- CO8756_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8756_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.160188
1985 2.999615
1986 2.990690
1987 3.483539
1988 3.821189
1989 3.695634
1990 3.177951
1991 3.774556
1992 2.954804
1993 2.786229
1994 3.132801
1995 2.943160
1996 3.420263
1997 3.773308
1998 2.210190
1999 4.035163
2000 3.458061
2001 3.343180
2002 3.140070
2003 3.572253
2004 3.422953
2005 3.577687
2006 3.287025
2007 3.748500
2008 4.291172
2009 4.260819
2010 4.029531
2011 3.548875
2012 3.863075
2013 4.991993
2014 4.468761
2015 3.181396
2016 3.759473
2017 4.348649
2018 2.815433
2019 3.898060
2020 4.884858
2021 3.814800
2022 4.678236 
ggplot(CO8756_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8756 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO8756_sd_mk_spring <- mk.test(CO8756_standard_dev_all_spring$sd_2)
print(CO8756_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO8756_standard_dev_all_spring$sd_2
## z = 2.9516, n = 39, p-value = 0.003161
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  245.0000000 6833.6666667    0.3306343
CO8756_sd_sens_spring <- sens.slope(CO8756_standard_dev_all_spring$sd_2)
print(CO8756_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO8756_standard_dev_all_spring$sd_2
## z = 2.9516, n = 39, p-value = 0.003161
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.009709191 0.041588410
## sample estimates:
## Sen's slope 
##  0.02775287

Fall temperature standard deviation

CO8756_standard_dev_all_fall <- CO8756_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO8756_standard_dev_all_fall <- CO8756_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO8756_standard_dev_all_fall <- CO8756_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8756_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.711644
1985 4.009774
1986 2.885793
1987 2.566292
1988 2.653328
1989 2.983588
1990 3.092729
1991 3.766456
1992 3.587660
1993 3.067939
1994 2.736613
1995 2.992300
1996 3.822336
1997 3.947173
1998 3.312247
1999 3.068790
2000 3.111605
2001 2.926958
2002 3.105300
2003 2.351483
2004 3.430668
2005 2.690507
2006 3.944444
2007 3.785637
2008 3.495235
2009 3.590275
2010 4.142315
2011 3.184582
2012 3.785243
2013 4.141586
2014 3.744568
2015 3.348521
2016 3.181286
2017 4.087696
2018 3.979366
2019 4.016837
2020 5.902965
2021 5.968943
2022 3.266773 
ggplot(CO8756_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8756 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO8756_sd_mk_fall <- mk.test(CO8756_standard_dev_all_fall$sd_2)
print(CO8756_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO8756_standard_dev_all_fall$sd_2
## z = 3.7016, n = 39, p-value = 0.0002142
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##           S        varS         tau 
##  307.000000 6833.666667    0.414305
CO8756_sd_sens_fall <- sens.slope(CO8756_standard_dev_all_fall$sd_2)
print(CO8756_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO8756_standard_dev_all_fall$sd_2
## z = 3.7016, n = 39, p-value = 0.0002142
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.01388079 0.04741831
## sample estimates:
## Sen's slope 
##  0.03186578

KREMMLING CO4664

Detrending Data

#average water year temperature

CO4664_yearly_wy_aver <- COOP_stations_clean %>% 
  filter(station == "CO4664") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO4664_daily_wy_aver <- CO4664_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

CO4664_daily_wy_aver <- CO4664_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO4664_daily_wy_aver$aver_day_temp))

#str(CO4664_daily_wy_aver)
# try to show all years as means. 
CO4664_daily_wy_aver2 <- CO4664_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO4664_daily_wy_aver2$date_temp <- signif(CO4664_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO4664_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO4664_standard_dev <- CO4664_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO4664_standard_dev$residual)
## [1] 2.341221e-16

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO4664_standard_dev_all <- CO4664_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO4664_standard_dev_all <- CO4664_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO4664_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.629316
1985 3.925765
1986 4.391539
1987 3.578252
1988 4.444457
1989 3.826456
1990 3.544877
1991 3.982667
1992 3.681005
1993 4.435733
1994 4.106463
1995 4.375472
1996 4.481437
1997 4.383461
1998 3.764676
1999 4.774303
2000 3.998953
2001 4.002248
2002 3.988324
2003 4.199619
2004 4.162573
2005 4.256445
2006 4.619417
2007 4.503859
2008 4.048740
2009 4.131050
2010 4.234699
2011 4.316928
2012 3.823873
2013 5.389521
2014 4.349641
2015 4.178394
2016 4.021211
2017 4.671163
2018 3.804800
2019 4.020858
2020 4.611534
2021 4.028081
2022 5.017778 
ggplot(CO4664_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO4664_sd_mk <- mk.test(CO4664_standard_dev_all$sd_2)
print(CO4664_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO4664_standard_dev_all$sd_2
## z = 1.3307, n = 39, p-value = 0.1833
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  111.0000000 6833.6666667    0.1497976
CO4664_sd_sens <- sens.slope(CO4664_standard_dev_all$sd_2)
print(CO4664_sd_sens)
## 
##  Sen's slope
## 
## data:  CO4664_standard_dev_all$sd_2
## z = 1.3307, n = 39, p-value = 0.1833
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.004123049  0.019440137
## sample estimates:
## Sen's slope 
## 0.007890612

Summer temperature standard deviation

CO4664_standard_dev_all_summer <- CO4664_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO4664_standard_dev_all_summer <- CO4664_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO4664_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.233705
1985 2.165484
1986 2.065486
1987 2.072465
1988 2.047536
1989 1.863139
1990 2.248649
1991 1.791980
1992 2.316931
1993 2.477050
1994 2.244732
1995 2.876224
1996 1.740202
1997 2.421018
1998 2.654645
1999 1.964120
2000 2.248942
2001 2.608910
2002 2.114431
2003 2.283702
2004 2.632135
2005 2.388002
2006 2.272657
2007 2.238510
2008 2.416289
2009 1.958795
2010 2.476500
2011 1.876893
2012 1.863630
2013 2.146713
2014 2.033889
2015 2.350329
2016 2.159097
2017 1.900767
2018 2.177626
2019 2.485504
2020 2.320842
2021 2.252179
2022 2.294532 
ggplot(CO4664_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO4664 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO4664_sd_mk_summer <- mk.test(CO4664_standard_dev_all_summer$sd_2)
print(CO4664_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO4664_standard_dev_all_summer$sd_2
## z = 0.60484, n = 39, p-value = 0.5453
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 5.100000e+01 6.833667e+03 6.882591e-02
CO4664_sd_sens_summer <- sens.slope(CO4664_standard_dev_all_summer$sd_2)
print(CO4664_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO4664_standard_dev_all_summer$sd_2
## z = 0.60484, n = 39, p-value = 0.5453
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.006626342  0.008864937
## sample estimates:
## Sen's slope 
## 0.002072301

Winter temperature standard deviation

CO4664_standard_dev_all_winter <- CO4664_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO4664_standard_dev_all_winter <- CO4664_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO4664_standard_dev_all_winter <- CO4664_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO4664_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 6.178206
1985 4.884856
1986 6.007381
1987 4.561534
1988 5.758177
1989 4.640323
1990 4.052267
1991 5.138970
1992 4.236067
1993 5.818980
1994 5.501827
1995 5.830597
1996 5.998059
1997 5.682426
1998 4.649196
1999 6.118840
2000 4.969854
2001 4.862641
2002 5.104401
2003 5.542611
2004 5.086410
2005 5.731239
2006 6.162130
2007 6.048120
2008 5.151940
2009 5.118672
2010 5.130619
2011 5.824734
2012 4.832756
2013 7.114688
2014 5.591827
2015 5.619818
2016 5.383425
2017 5.942280
2018 4.646987
2019 5.029837
2020 5.543155
2021 4.427400
2022 6.736750 
ggplot(CO4664_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO4664 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO4664_sd_mk_winter <- mk.test(CO4664_standard_dev_all_winter$sd_2)
print(CO4664_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO4664_standard_dev_all_winter$sd_2
## z = 0.55646, n = 39, p-value = 0.5779
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   47.0000000 6833.6666667    0.0634278
CO4664_sd_sens_winter <- sens.slope(CO4664_standard_dev_all_winter$sd_2)
print(CO4664_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO4664_standard_dev_all_winter$sd_2
## z = 0.55646, n = 39, p-value = 0.5779
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01336425  0.02973909
## sample estimates:
## Sen's slope 
##  0.00618425

Spring and Fall temperature standard deviation

CO4664_standard_dev_all_spring <- CO4664_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO4664_standard_dev_all_spring <- CO4664_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO4664_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.673611
1985 2.931194
1986 3.031972
1987 2.794903
1988 3.359923
1989 3.576750
1990 2.809345
1991 3.555546
1992 2.749745
1993 2.877329
1994 3.064787
1995 2.506119
1996 3.516907
1997 3.519984
1998 2.628799
1999 4.045450
2000 3.765940
2001 3.577176
2002 3.275739
2003 3.651328
2004 3.225955
2005 3.124731
2006 3.148553
2007 3.551324
2008 3.336745
2009 3.644830
2010 3.763705
2011 3.594960
2012 3.525322
2013 4.306641
2014 4.351558
2015 2.628487
2016 3.068566
2017 4.106248
2018 2.598629
2019 3.589900
2020 3.892748
2021 3.548363
2022 3.911897 
ggplot(CO4664_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO4664 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO4664_sd_mk_spring <- mk.test(CO4664_standard_dev_all_spring$sd_2)
print(CO4664_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO4664_standard_dev_all_spring$sd_2
## z = 2.2984, n = 39, p-value = 0.02154
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  191.0000000 6833.6666667    0.2577598
CO4664_sd_sens_spring <- sens.slope(CO4664_standard_dev_all_spring$sd_2)
print(CO4664_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO4664_standard_dev_all_spring$sd_2
## z = 2.2984, n = 39, p-value = 0.02154
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.002162987 0.030412117
## sample estimates:
## Sen's slope 
##  0.01690692

Fall temperature standard deviation

CO4664_standard_dev_all_fall <- CO4664_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO4664_standard_dev_all_fall <- CO4664_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO4664_standard_dev_all_fall <- CO4664_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO4664_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.694545
1985 3.348074
1986 2.841655
1987 2.434565
1988 3.147555
1989 2.842991
1990 2.922646
1991 3.290714
1992 3.708803
1993 3.466340
1994 2.778215
1995 3.102763
1996 3.688956
1997 3.530037
1998 3.518928
1999 3.194830
2000 3.082765
2001 2.775091
2002 2.871231
2003 2.518007
2004 3.631037
2005 2.560336
2006 3.222795
2007 3.099360
2008 2.943519
2009 3.280135
2010 3.702819
2011 2.662267
2012 3.364123
2013 3.375674
2014 3.343279
2015 2.621939
2016 2.417112
2017 3.312480
2018 3.110671
2019 3.253837
2020 4.807376
2021 5.284076
2022 2.758460 
ggplot(CO4664_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO4664 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO4664_sd_mk_fall <- mk.test(CO4664_standard_dev_all_fall$sd_2)
print(CO4664_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO4664_standard_dev_all_fall$sd_2
## z = 0.79839, n = 39, p-value = 0.4246
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 6.700000e+01 6.833667e+03 9.041835e-02
CO4664_sd_sens_fall <- sens.slope(CO4664_standard_dev_all_fall$sd_2)
print(CO4664_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO4664_standard_dev_all_fall$sd_2
## z = 0.79839, n = 39, p-value = 0.4246
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.008501993  0.019778494
## sample estimates:
## Sen's slope 
## 0.006715174

CRAIG 4 SW CO1932

Detrending Data

#average water year temperature

CO1932_yearly_wy_aver <- COOP_stations_clean %>% 
  filter(station == "CO1932") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO1932_daily_wy_aver <- CO1932_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

CO1932_daily_wy_aver <- CO1932_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO1932_daily_wy_aver$aver_day_temp))

#str(CO1932_daily_wy_aver)
# try to show all years as means. 
CO1932_daily_wy_aver2 <- CO1932_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO1932_daily_wy_aver2$date_temp <- signif(CO1932_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO1932_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO1932_standard_dev <- CO1932_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO1932_standard_dev$residual)
## [1] 1.452822e-16

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO1932_standard_dev_all <- CO1932_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO1932_standard_dev_all <- CO1932_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO1932_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.233943
1985 4.492398
1986 4.288049
1987 4.086941
1988 4.110339
1989 4.952326
1990 3.781078
1991 4.332661
1992 4.274896
1993 3.975801
1994 3.765651
1995 4.187504
1996 3.988335
1997 4.314615
1998 3.695730
1999 4.074580
2000 3.881989
2001 3.805798
2002 4.114684
2003 3.814501
2004 4.064477
2005 4.045794
2006 4.361321
2007 4.396260
2008 4.067821
2009 4.148666
2010 4.194158
2011 4.462111
2012 3.624319
2013 4.797391
2014 4.293946
2015 4.598871
2016 4.014741
2017 4.484346
2018 4.031234
2019 4.108593
2020 4.610221
2021 4.224450
2022 4.470230 
ggplot(CO1932_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO1932_sd_mk <- mk.test(CO1932_standard_dev_all$sd_2)
print(CO1932_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO1932_standard_dev_all$sd_2
## z = 1.0887, n = 39, p-value = 0.2763
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##           S        varS         tau 
##   91.000000 6833.666667    0.122807
CO1932_sd_sens <- sens.slope(CO1932_standard_dev_all$sd_2)
print(CO1932_sd_sens)
## 
##  Sen's slope
## 
## data:  CO1932_standard_dev_all$sd_2
## z = 1.0887, n = 39, p-value = 0.2763
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.00360703  0.01394344
## sample estimates:
## Sen's slope 
## 0.005228625

Summer temperature standard deviation

CO1932_standard_dev_all_summer <- CO1932_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO1932_standard_dev_all_summer <- CO1932_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO1932_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.645536
1985 2.790609
1986 2.302503
1987 2.598101
1988 2.362997
1989 2.485319
1990 2.729863
1991 2.001142
1992 2.732537
1993 2.705879
1994 1.973762
1995 2.705922
1996 1.766488
1997 2.362222
1998 2.806903
1999 1.784271
2000 2.240505
2001 2.463228
2002 2.564459
2003 2.360163
2004 2.627324
2005 2.703909
2006 2.351098
2007 2.370944
2008 2.860657
2009 2.175469
2010 2.491530
2011 2.011191
2012 2.292131
2013 1.958284
2014 2.498162
2015 2.510988
2016 2.189942
2017 2.104386
2018 2.546445
2019 2.492025
2020 2.675709
2021 2.778538
2022 2.356130 
ggplot(CO1932_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO1932 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO1932_sd_mk_summer <- mk.test(CO1932_standard_dev_all_summer$sd_2)
print(CO1932_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO1932_standard_dev_all_summer$sd_2
## z = -0.62904, n = 39, p-value = 0.5293
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -53.00000000 6833.66666667   -0.07152497
CO1932_sd_sens_summer <- sens.slope(CO1932_standard_dev_all_summer$sd_2)
print(CO1932_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO1932_standard_dev_all_summer$sd_2
## z = -0.62904, n = 39, p-value = 0.5293
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01059128  0.00608767
## sample estimates:
## Sen's slope 
## -0.00202955

Winter temperature standard deviation

CO1932_standard_dev_all_winter <- CO1932_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO1932_standard_dev_all_winter <- CO1932_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO1932_standard_dev_all_winter <- CO1932_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO1932_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.747174
1985 4.991136
1986 5.276625
1987 4.964778
1988 4.927283
1989 6.305333
1990 4.373717
1991 5.653904
1992 4.751119
1993 4.812670
1994 4.713050
1995 5.251216
1996 4.721374
1997 4.977345
1998 4.456422
1999 4.683508
2000 4.328575
2001 4.367117
2002 5.091952
2003 4.447521
2004 4.779659
2005 4.892874
2006 5.482335
2007 5.671181
2008 4.639198
2009 5.081722
2010 4.958611
2011 5.878069
2012 4.338114
2013 5.837190
2014 5.197445
2015 6.094602
2016 4.950442
2017 5.592113
2018 4.817253
2019 4.845713
2020 4.489499
2021 4.263175
2022 5.673973 
ggplot(CO1932_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO1932 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO1932_sd_mk_winter <- mk.test(CO1932_standard_dev_all_winter$sd_2)
print(CO1932_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO1932_standard_dev_all_winter$sd_2
## z = 0.24194, n = 39, p-value = 0.8088
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 2.100000e+01 6.833667e+03 2.834008e-02
CO1932_sd_sens_winter <- sens.slope(CO1932_standard_dev_all_winter$sd_2)
print(CO1932_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO1932_standard_dev_all_winter$sd_2
## z = 0.24194, n = 39, p-value = 0.8088
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.0125442  0.0187884
## sample estimates:
## Sen's slope 
##  0.00150348

Spring and Fall temperature standard deviation

CO1932_standard_dev_all_spring <- CO1932_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO1932_standard_dev_all_spring <- CO1932_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO1932_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.755314
1985 3.903831
1986 3.767537
1987 3.786244
1988 4.149960
1989 4.382411
1990 3.541878
1991 3.737712
1992 3.246495
1993 3.353646
1994 3.697772
1995 2.805358
1996 3.765897
1997 4.417088
1998 2.615912
1999 3.871307
2000 3.962378
2001 3.673094
2002 3.227299
2003 4.313896
2004 3.687416
2005 3.373838
2006 3.334905
2007 3.633898
2008 4.100934
2009 4.073092
2010 4.136656
2011 3.740589
2012 3.686838
2013 5.182217
2014 4.274755
2015 3.142451
2016 3.515847
2017 3.903040
2018 3.302935
2019 3.803930
2020 4.489189
2021 4.037137
2022 4.193143 
ggplot(CO1932_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO1932 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO1932_sd_mk_spring <- mk.test(CO1932_standard_dev_all_spring$sd_2)
print(CO1932_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO1932_standard_dev_all_spring$sd_2
## z = 0.31452, n = 39, p-value = 0.7531
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 2.700000e+01 6.833667e+03 3.643725e-02
CO1932_sd_sens_spring <- sens.slope(CO1932_standard_dev_all_spring$sd_2)
print(CO1932_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO1932_standard_dev_all_spring$sd_2
## z = 0.31452, n = 39, p-value = 0.7531
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01089316  0.01809421
## sample estimates:
## Sen's slope 
## 0.003559963

Fall temperature standard deviation

CO1932_standard_dev_all_fall <- CO1932_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO1932_standard_dev_all_fall <- CO1932_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO1932_standard_dev_all_fall <- CO1932_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO1932_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.564447
1985 4.544833
1986 3.932172
1987 3.335138
1988 3.305449
1989 3.413709
1990 3.681264
1991 3.637236
1992 4.537751
1993 3.466026
1994 2.960881
1995 2.974595
1996 3.965793
1997 4.378966
1998 3.330774
1999 3.313114
2000 3.748205
2001 3.064185
2002 3.338009
2003 2.775139
2004 4.015253
2005 2.644792
2006 4.228341
2007 3.355968
2008 3.697894
2009 3.264461
2010 4.113869
2011 2.992072
2012 3.067413
2013 3.925543
2014 4.047597
2015 3.106753
2016 3.290819
2017 4.150447
2018 3.991989
2019 3.983640
2020 6.218528
2021 5.584914
2022 3.271232 
ggplot(CO1932_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO1932 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO1932_sd_mk_fall <- mk.test(CO1932_standard_dev_all_fall$sd_2)
print(CO1932_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO1932_standard_dev_all_fall$sd_2
## z = 0.53226, n = 39, p-value = 0.5945
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 4.500000e+01 6.833667e+03 6.072874e-02
CO1932_sd_sens_fall <- sens.slope(CO1932_standard_dev_all_fall$sd_2)
print(CO1932_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO1932_standard_dev_all_fall$sd_2
## z = 0.53226, n = 39, p-value = 0.5945
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01314052  0.02480373
## sample estimates:
## Sen's slope 
## 0.003765486

STEAMBOAT SPRINGS CO7936

Detrending Data

#average water year temperature

CO7936_yearly_wy_aver <- COOP_stations_clean %>%
  filter(station == "CO7936") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO7936_daily_wy_aver <- CO7936_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

CO7936_daily_wy_aver <- CO7936_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO7936_daily_wy_aver$aver_day_temp))

#str(CO7936_daily_wy_aver)
# try to show all years as means. 
CO7936_daily_wy_aver2 <- CO7936_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO7936_daily_wy_aver2$date_temp <- signif(CO7936_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO7936_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO7936_standard_dev <- CO7936_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO7936_standard_dev$residual)
## [1] 2.347603e-17

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO7936_standard_dev_all <- CO7936_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO7936_standard_dev_all <- CO7936_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7936_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.839061
1985 3.929760
1986 3.940912
1987 3.447446
1988 3.906482
1989 4.292675
1990 3.620956
1991 4.174814
1992 3.963401
1993 4.002728
1994 3.707582
1995 3.970729
1996 3.932423
1997 4.046533
1998 3.678931
1999 4.159525
2000 3.665091
2001 4.138010
2002 4.172440
2003 3.723505
2004 3.878909
2005 3.683221
2006 4.118251
2007 4.041718
2008 3.717350
2009 3.745580
2010 3.935927
2011 4.188137
2012 3.315756
2013 4.330880
2014 3.861902
2015 4.093439
2016 3.564639
2017 4.395071
2018 3.647586
2019 3.969781
2020 4.359443
2021 3.804308
2022 4.155187 
#CALLING THIS something different
CO7936_all_V2 <- ggplot(CO7936_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')


print(CO7936_all_V2)

CO7936_sd_mk <- mk.test(CO7936_standard_dev_all$sd_2)
print(CO7936_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO7936_standard_dev_all$sd_2
## z = 0.87097, n = 39, p-value = 0.3838
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 7.300000e+01 6.833667e+03 9.851552e-02
CO7936_sd_sens <- sens.slope(CO7936_standard_dev_all$sd_2)
print(CO7936_sd_sens)
## 
##  Sen's slope
## 
## data:  CO7936_standard_dev_all$sd_2
## z = 0.87097, n = 39, p-value = 0.3838
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.004751728  0.011654360
## sample estimates:
## Sen's slope 
## 0.003841881

Summer temperature standard deviation

CO7936_standard_dev_all_summer <- CO7936_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO7936_standard_dev_all_summer <- CO7936_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7936_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 1.962215
1985 2.126840
1986 1.814465
1987 1.950477
1988 2.004603
1989 2.223426
1990 2.568322
1991 1.747480
1992 2.340906
1993 2.554491
1994 2.181100
1995 2.758238
1996 1.772840
1997 2.350319
1998 2.386414
1999 1.831799
2000 2.222008
2001 2.646216
2002 2.533728
2003 2.402936
2004 2.497338
2005 2.235540
2006 2.136381
2007 2.366067
2008 2.562213
2009 1.938279
2010 2.134068
2011 1.848619
2012 1.993840
2013 1.931257
2014 2.224286
2015 2.314251
2016 2.007504
2017 1.886006
2018 2.283866
2019 2.386876
2020 2.530220
2021 2.461403
2022 2.262040 
ggplot(CO7936_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7936 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO7936_sd_mk_summer <- mk.test(CO7936_standard_dev_all_summer$sd_2)
print(CO7936_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO7936_standard_dev_all_summer$sd_2
## z = 0.84678, n = 39, p-value = 0.3971
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 7.100000e+01 6.833667e+03 9.581646e-02
CO7936_sd_sens_summer <- sens.slope(CO7936_standard_dev_all_summer$sd_2)
print(CO7936_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO7936_standard_dev_all_summer$sd_2
## z = 0.84678, n = 39, p-value = 0.3971
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.005127418  0.011525156
## sample estimates:
## Sen's slope 
## 0.003416092

Winter temperature standard deviation

CO7936_standard_dev_all_winter <- CO7936_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO7936_standard_dev_all_winter <- CO7936_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO7936_standard_dev_all_winter <- CO7936_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7936_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.892874
1985 5.034645
1986 5.246769
1987 4.209801
1988 4.943908
1989 5.585204
1990 4.300953
1991 5.456982
1992 4.579049
1993 5.046461
1994 4.683210
1995 4.953747
1996 4.816397
1997 4.691776
1998 4.646675
1999 5.192548
2000 4.178326
2001 4.871604
2002 5.069140
2003 4.574573
2004 4.705222
2005 4.589251
2006 5.142351
2007 5.041724
2008 4.251188
2009 4.523145
2010 4.693132
2011 5.548294
2012 3.971197
2013 5.367923
2014 4.640035
2015 5.130570
2016 4.250227
2017 5.457668
2018 4.384247
2019 4.671086
2020 4.587241
2021 4.047675
2022 5.184049 
ggplot(CO7936_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7936 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO7936_sd_mk_winter <- mk.test(CO7936_standard_dev_all_winter$sd_2)
print(CO7936_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO7936_standard_dev_all_winter$sd_2
## z = -1.0645, n = 39, p-value = 0.2871
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##           S        varS         tau 
##  -89.000000 6833.666667   -0.120108
CO7936_sd_sens_winter <- sens.slope(CO7936_standard_dev_all_winter$sd_2)
print(CO7936_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO7936_standard_dev_all_winter$sd_2
## z = -1.0645, n = 39, p-value = 0.2871
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.02008272  0.00584542
## sample estimates:
## Sen's slope 
##  -0.0082035

Spring and Fall temperature standard deviation

CO7936_standard_dev_all_spring <- CO7936_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO7936_standard_dev_all_spring <- CO7936_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7936_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.790439
1985 2.902924
1986 3.083814
1987 3.127936
1988 3.522284
1989 3.303794
1990 3.599343
1991 3.879164
1992 3.005685
1993 3.409430
1994 3.270398
1995 2.883955
1996 3.577860
1997 3.936294
1998 2.590320
1999 4.259020
2000 3.872380
2001 3.209207
2002 3.135469
2003 3.878654
2004 3.394535
2005 3.488069
2006 3.188137
2007 3.682602
2008 3.845791
2009 3.726660
2010 3.891196
2011 3.164650
2012 3.515949
2013 4.451575
2014 3.856463
2015 2.999571
2016 3.288434
2017 4.052723
2018 3.108385
2019 3.745441
2020 4.197006
2021 3.729601
2022 4.065043 
ggplot(CO7936_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7936 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO7936_sd_mk_spring <- mk.test(CO7936_standard_dev_all_spring$sd_2)
print(CO7936_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO7936_standard_dev_all_spring$sd_2
## z = 2.129, n = 39, p-value = 0.03325
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  177.0000000 6833.6666667    0.2388664
CO7936_sd_sens_spring <- sens.slope(CO7936_standard_dev_all_spring$sd_2)
print(CO7936_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO7936_standard_dev_all_spring$sd_2
## z = 2.129, n = 39, p-value = 0.03325
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.001778762 0.026982490
## sample estimates:
## Sen's slope 
##  0.01375395

Fall temperature standard deviation

CO7936_standard_dev_all_fall <- CO7936_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO7936_standard_dev_all_fall <- CO7936_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO7936_standard_dev_all_fall <- CO7936_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7936_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.833969
1985 3.402018
1986 2.846760
1987 2.618908
1988 2.675820
1989 2.778285
1990 2.904594
1991 3.353640
1992 3.826229
1993 3.199968
1994 2.755089
1995 3.064669
1996 3.902089
1997 4.425399
1998 3.173550
1999 3.439840
2000 3.339817
2001 2.842223
2002 3.127491
2003 2.445499
2004 3.560674
2005 2.400901
2006 3.742896
2007 3.233108
2008 3.277337
2009 2.915478
2010 3.936816
2011 3.199677
2012 2.824785
2013 3.620541
2014 3.698383
2015 3.091351
2016 3.083811
2017 3.705909
2018 3.410984
2019 3.984335
2020 5.413457
2021 4.738721
2022 3.105386 
ggplot(CO7936_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7936 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO7936_sd_mk_fall <- mk.test(CO7936_standard_dev_all_fall$sd_2)
print(CO7936_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO7936_standard_dev_all_fall$sd_2
## z = 2.5403, n = 39, p-value = 0.01107
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  211.0000000 6833.6666667    0.2847503
CO7936_sd_sens_fall <- sens.slope(CO7936_standard_dev_all_fall$sd_2)
print(CO7936_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO7936_standard_dev_all_fall$sd_2
## z = 2.5403, n = 39, p-value = 0.01107
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.00585575 0.03623337
## sample estimates:
## Sen's slope 
##  0.01847364

Southern Stations

HERMIT 7 ESE CO3951

Detrending Data

#average water year temperature

CO3951_yearly_wy_aver <- COOP_stations_clean %>% 
  filter(station == "CO3951") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO3951_daily_wy_aver <- CO3951_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

CO3951_daily_wy_aver <- CO3951_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO3951_daily_wy_aver$aver_day_temp))

#str(CO3951_daily_wy_aver)
# try to show all years as means. 
CO3951_daily_wy_aver2 <- CO3951_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO3951_daily_wy_aver2$date_temp <- signif(CO3951_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO3951_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO3951_standard_dev <- CO3951_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO3951_standard_dev$residual)
## [1] -2.03745e-16

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO3951_standard_dev_all <- CO3951_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO3951_standard_dev_all <- CO3951_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3951_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.653442
1985 3.896927
1986 2.865464
1987 3.240336
1988 3.240178
1989 3.950884
1990 3.319987
1991 3.981485
1992 3.275238
1993 3.848168
1994 3.671410
1995 3.526380
1996 3.554982
1997 3.335595
1998 3.115427
1999 3.853665
2000 3.469245
2001 3.574442
2002 3.357829
2003 3.384614
2004 4.035331
2005 3.434813
2006 3.654846
2007 3.634470
2008 4.205183
2009 3.578311
2010 3.433010
2011 3.919643
2012 3.485099
2013 3.944313
2014 3.535969
2015 3.615229
2016 3.743967
2017 4.003013
2018 3.285469
2019 3.113867
2020 3.118208
2021 2.781827
2022 2.949074 
ggplot(CO3951_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3951_sd_mk <- mk.test(CO3951_standard_dev_all$sd_2)
print(CO3951_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO3951_standard_dev_all$sd_2
## z = -0.33871, n = 39, p-value = 0.7348
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  -29.0000000 6833.6666667   -0.0391363
CO3951_sd_sens <- sens.slope(CO3951_standard_dev_all$sd_2)
print(CO3951_sd_sens)
## 
##  Sen's slope
## 
## data:  CO3951_standard_dev_all$sd_2
## z = -0.33871, n = 39, p-value = 0.7348
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.014902241  0.009629078
## sample estimates:
##  Sen's slope 
## -0.002428287

Summer temperature standard deviation

CO3951_standard_dev_all_summer <- CO3951_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO3951_standard_dev_all_summer <- CO3951_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3951_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 1.611452
1985 1.821910
1986 1.678949
1987 1.933003
1988 2.554016
1989 1.992681
1990 2.007483
1991 1.495795
1992 1.761612
1993 2.265070
1994 1.538173
1995 2.718106
1996 1.883564
1997 1.927223
1998 2.084648
1999 1.905750
2000 1.375707
2001 2.026183
2002 1.803919
2003 2.098167
2004 2.937211
2005 2.205640
2006 2.017768
2007 1.786142
2008 2.048327
2009 2.367984
2010 2.129345
2011 1.740053
2012 1.630538
2013 1.822416
2014 1.755879
2015 2.267538
2016 2.287124
2017 2.241756
2018 1.205778
2019 1.547015
2020 1.656319
2021 1.733652
2022 1.517842 
ggplot(CO3951_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3951 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO3951_sd_mk_summer <- mk.test(CO3951_standard_dev_all_summer$sd_2)
print(CO3951_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO3951_standard_dev_all_summer$sd_2
## z = -0.3871, n = 39, p-value = 0.6987
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -33.00000000 6833.66666667   -0.04453441
CO3951_sd_sens_summer <- sens.slope(CO3951_standard_dev_all_summer$sd_2)
print(CO3951_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO3951_standard_dev_all_summer$sd_2
## z = -0.3871, n = 39, p-value = 0.6987
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01273478  0.00984424
## sample estimates:
##  Sen's slope 
## -0.002463427

Winter temperature standard deviation

CO3951_standard_dev_all_winter <- CO3951_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO3951_standard_dev_all_winter <- CO3951_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO3951_standard_dev_all_winter <- CO3951_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3951_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.376898
1985 5.110075
1986 3.498172
1987 3.818331
1988 4.041878
1989 4.938470
1990 4.241458
1991 4.221924
1992 3.429547
1993 3.486827
1994 4.660710
1995 4.224296
1996 4.275101
1997 4.113673
1998 3.833033
1999 3.655830
2000 4.338066
2001 4.570178
2002 4.265218
2003 4.059037
2004 4.466516
2005 4.331086
2006 4.668315
2007 4.638112
2008 4.782396
2009 4.210479
2010 3.752505
2011 5.130373
2012 4.319224
2013 5.227346
2014 4.127587
2015 4.363078
2016 4.839900
2017 4.302162
2018 3.970735
2019 3.971246
2020 3.748711
2021 3.585847
2022 3.568386 
ggplot(CO3951_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3951 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO3951_sd_mk_winter <- mk.test(CO3951_standard_dev_all_winter$sd_2)
print(CO3951_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO3951_standard_dev_all_winter$sd_2
## z = -0.024194, n = 39, p-value = 0.9807
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
## -3.000000e+00  6.833667e+03 -4.048583e-03
CO3951_sd_sens_winter <- sens.slope(CO3951_standard_dev_all_winter$sd_2)
print(CO3951_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO3951_standard_dev_all_winter$sd_2
## z = -0.024194, n = 39, p-value = 0.9807
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01642490  0.01504415
## sample estimates:
##   Sen's slope 
## -0.0004458205

Spring and Fall temperature standard deviation

CO3951_standard_dev_all_spring <- CO3951_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO3951_standard_dev_all_spring <- CO3951_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3951_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.974859
1985 2.266381
1986 2.596851
1987 3.001591
1988 2.445149
1989 2.767623
1990 1.886559
1991 3.959482
1992 2.160393
1993 3.658601
1994 2.964015
1995 3.233875
1996 2.412861
1997 3.061968
1998 2.690651
1999 3.737786
2000 3.433883
2001 3.374779
2002 2.899746
2003 3.299711
2004 2.881078
2005 3.221468
2006 2.705256
2007 3.684370
2008 3.677992
2009 2.964750
2010 3.771933
2011 3.799636
2012 3.457868
2013 3.205095
2014 3.783368
2015 3.229056
2016 3.279492
2017 3.376814
2018 2.461386
2019 2.998422
2020 2.402555
2021 2.304718
2022 2.681835 
ggplot(CO3951_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3951 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO3951_sd_mk_spring <- mk.test(CO3951_standard_dev_all_spring$sd_2)
print(CO3951_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO3951_standard_dev_all_spring$sd_2
## z = 0.55646, n = 39, p-value = 0.5779
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   47.0000000 6833.6666667    0.0634278
CO3951_sd_sens_spring <- sens.slope(CO3951_standard_dev_all_spring$sd_2)
print(CO3951_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO3951_standard_dev_all_spring$sd_2
## z = 0.55646, n = 39, p-value = 0.5779
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01365519  0.02743630
## sample estimates:
## Sen's slope 
##  0.00580447

Fall temperature standard deviation

CO3951_standard_dev_all_fall <- CO3951_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO3951_standard_dev_all_fall <- CO3951_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO3951_standard_dev_all_fall <- CO3951_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3951_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.172938
1985 3.163365
1986 2.419472
1987 2.464200
1988 2.291031
1989 2.178750
1990 2.522609
1991 2.014982
1992 2.636411
1993 2.510656
1994 3.617128
1995 2.567791
1996 2.720100
1997 3.006096
1998 2.794716
1999 3.604555
2000 2.806227
2001 2.264470
2002 2.361697
2003 2.772334
2004 4.259083
2005 2.367137
2006 2.913526
2007 2.555338
2008 3.046940
2009 2.974680
2010 3.634692
2011 2.705310
2012 2.608423
2013 2.743032
2014 3.559448
2015 2.234117
2016 2.519028
2017 2.880449
2018 1.745936
2019 2.468167
2020 3.075222
2021 2.173653
2022 2.402625 
ggplot(CO3951_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3951 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO3951_sd_mk_fall <- mk.test(CO3951_standard_dev_all_fall$sd_2)
print(CO3951_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO3951_standard_dev_all_fall$sd_2
## z = 0.70162, n = 39, p-value = 0.4829
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 5.900000e+01 6.833667e+03 7.962213e-02
CO3951_sd_sens_fall <- sens.slope(CO3951_standard_dev_all_fall$sd_2)
print(CO3951_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO3951_standard_dev_all_fall$sd_2
## z = 0.70162, n = 39, p-value = 0.4829
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.009813804  0.019286754
## sample estimates:
## Sen's slope 
## 0.005077217

PAGOSA SPRINGS CO6258

Detrending Data

#average water year temperature

CO6258_yearly_wy_aver <- COOP_stations_clean %>% 
  filter(station == "CO6258") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO6258_daily_wy_aver <- CO6258_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

CO6258_daily_wy_aver <- CO6258_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO6258_daily_wy_aver$aver_day_temp))

#str(CO6258_daily_wy_aver)
# try to show all years as means. 
CO6258_daily_wy_aver2 <- CO6258_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO6258_daily_wy_aver2$date_temp <- signif(CO6258_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO6258_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO6258_standard_dev <- CO6258_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO6258_standard_dev$residual)
## [1] 3.510867e-16

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO6258_standard_dev_all <- CO6258_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO6258_standard_dev_all <- CO6258_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO6258_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.318666
1985 3.548810
1986 3.585376
1987 3.241446
1988 3.268235
1989 3.405561
1990 3.140864
1991 3.359892
1992 3.211218
1993 3.204075
1994 2.753097
1995 3.166953
1996 3.005127
1997 3.156048
1998 3.024570
1999 3.061679
2000 2.835479
2001 2.879689
2002 2.845815
2003 2.722025
2004 3.063822
2005 2.806196
2006 3.028786
2007 3.046606
2008 3.101611
2009 2.824292
2010 2.910644
2011 3.337556
2012 2.906860
2013 3.730961
2014 2.950127
2015 3.103518
2016 3.073977
2017 3.251814
2018 2.819956
2019 3.206256
2020 3.660927 
ggplot(CO6258_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO6258_sd_mk <- mk.test(CO6258_standard_dev_all$sd_2)
print(CO6258_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO6258_standard_dev_all$sd_2
## z = -1.7395, n = 37, p-value = 0.08195
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -134.0000000 5846.0000000   -0.2012012
CO6258_sd_sens <- sens.slope(CO6258_standard_dev_all$sd_2)
print(CO6258_sd_sens)
## 
##  Sen's slope
## 
## data:  CO6258_standard_dev_all$sd_2
## z = -1.7395, n = 37, p-value = 0.08195
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.014942634  0.001058464
## sample estimates:
##  Sen's slope 
## -0.008921116

Summer temperature standard deviation

CO6258_standard_dev_all_summer <- CO6258_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO6258_standard_dev_all_summer <- CO6258_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO6258_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 1.858786
1985 2.308794
1986 2.222994
1987 2.269465
1988 2.257263
1989 1.987246
1990 2.293365
1991 1.891664
1992 2.180190
1993 1.952882
1994 1.292074
1995 2.439176
1996 1.552245
1997 1.800557
1998 2.358568
1999 1.397148
2000 1.352726
2001 1.712864
2002 1.577063
2003 1.780712
2004 1.914234
2005 1.657315
2006 1.757491
2007 1.457080
2008 1.533523
2009 1.600800
2010 1.841829
2011 1.403970
2012 1.426995
2013 1.932821
2014 1.610427
2015 1.851429
2016 2.198329
2017 1.678228
2018 1.546923
2019 1.819680 
ggplot(CO6258_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO6258 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO6258_sd_mk_summer <- mk.test(CO6258_standard_dev_all_summer$sd_2)
print(CO6258_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO6258_standard_dev_all_summer$sd_2
## z = -2.4381, n = 36, p-value = 0.01476
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -180.0000000 5390.0000000   -0.2857143
CO6258_sd_sens_summer <- sens.slope(CO6258_standard_dev_all_summer$sd_2)
print(CO6258_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO6258_standard_dev_all_summer$sd_2
## z = -2.4381, n = 36, p-value = 0.01476
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.024358718 -0.003307517
## sample estimates:
## Sen's slope 
## -0.01496829

Winter temperature standard deviation

CO6258_standard_dev_all_winter <- CO6258_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO6258_standard_dev_all_winter <- CO6258_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO6258_standard_dev_all_winter <- CO6258_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO6258_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.024948
1985 4.137172
1986 4.246836
1987 3.803001
1988 3.913663
1989 4.296399
1990 3.869694
1991 4.245098
1992 3.698392
1993 3.904375
1994 3.360276
1995 3.394068
1996 3.576904
1997 3.740004
1998 3.418887
1999 3.277999
2000 3.566933
2001 3.473034
2002 3.522398
2003 3.159780
2004 3.754839
2005 3.499759
2006 3.698010
2007 3.912532
2008 3.945895
2009 3.381498
2010 3.024755
2011 4.358185
2012 3.565602
2013 4.671089
2014 3.245161
2015 3.717139
2016 3.697368
2017 3.912919
2018 3.437265
2019 3.689800
2020 3.343414 
ggplot(CO6258_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO6258 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO6258_sd_mk_winter <- mk.test(CO6258_standard_dev_all_winter$sd_2)
print(CO6258_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO6258_standard_dev_all_winter$sd_2
## z = -2.0272, n = 37, p-value = 0.04264
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -156.0000000 5846.0000000   -0.2342342
CO6258_sd_sens_winter <- sens.slope(CO6258_standard_dev_all_winter$sd_2)
print(CO6258_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO6258_standard_dev_all_winter$sd_2
## z = -2.0272, n = 37, p-value = 0.04264
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -2.218634e-02 -4.267337e-05
## sample estimates:
## Sen's slope 
## -0.01171235

Spring and Fall temperature standard deviation

CO6258_standard_dev_all_spring <- CO6258_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO6258_standard_dev_all_spring <- CO6258_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO6258_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.599752
1985 3.092177
1986 3.220708
1987 3.038607
1988 3.360173
1989 3.072899
1990 2.312228
1991 3.163667
1992 2.140066
1993 2.433322
1994 2.865133
1995 2.635910
1996 2.810729
1997 2.887475
1998 2.284744
1999 2.975571
2000 2.778467
2001 2.889011
2002 2.419112
2003 2.974263
2004 2.382293
2005 2.635896
2006 2.337629
2007 2.974535
2008 2.737540
2009 2.850171
2010 3.034482
2011 3.108309
2012 3.054762
2013 3.389964
2014 3.624183
2015 3.082779
2016 2.767349
2017 3.279508
2018 2.832477
2019 3.627504 
ggplot(CO6258_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO6258 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO6258_sd_mk_spring <- mk.test(CO6258_standard_dev_all_spring$sd_2)
print(CO6258_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO6258_standard_dev_all_spring$sd_2
## z = 0.74915, n = 36, p-value = 0.4538
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 5.600000e+01 5.390000e+03 8.888889e-02
CO6258_sd_sens_spring <- sens.slope(CO6258_standard_dev_all_spring$sd_2)
print(CO6258_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO6258_standard_dev_all_spring$sd_2
## z = 0.74915, n = 36, p-value = 0.4538
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.007360982  0.020287448
## sample estimates:
## Sen's slope 
## 0.007039707

Fall temperature standard deviation

CO6258_standard_dev_all_fall <- CO6258_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO6258_standard_dev_all_fall <- CO6258_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO6258_standard_dev_all_fall <- CO6258_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO6258_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.345071
1985 3.234957
1986 3.113612
1987 2.823486
1988 2.521678
1989 2.283309
1990 2.795173
1991 2.522469
1992 3.039080
1993 3.204105
1994 2.244921
1995 2.222479
1996 2.580649
1997 3.052734
1998 3.338668
1999 2.097198
2000 2.436414
2001 2.135227
2002 2.095117
2003 2.024071
2004 3.006352
2005 2.118445
2006 2.827557
2007 2.265079
2008 2.361061
2009 2.476929
2010 3.353483
2011 2.221133
2012 2.224251
2013 2.642040
2014 3.030271
2015 2.221308
2016 2.356646
2017 2.752698
2018 2.551065
2019 2.584246
2020 4.258011 
ggplot(CO6258_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO6258 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO6258_sd_mk_fall <- mk.test(CO6258_standard_dev_all_fall$sd_2)
print(CO6258_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO6258_standard_dev_all_fall$sd_2
## z = -0.4316, n = 37, p-value = 0.666
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -34.00000000 5846.00000000   -0.05105105
CO6258_sd_sens_fall <- sens.slope(CO6258_standard_dev_all_fall$sd_2)
print(CO6258_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO6258_standard_dev_all_fall$sd_2
## z = -0.4316, n = 37, p-value = 0.666
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01684571  0.01032787
## sample estimates:
##  Sen's slope 
## -0.002548867

SILVERTON CO7656

Detrending Data

#average water year temperature

CO7656_yearly_wy_aver <- COOP_stations_clean %>% 
  filter(station == "CO7656") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))

CO7656_temp_xts <- xts(CO7656_yearly_wy_aver$avg_T_c, order.by = CO7656_yearly_wy_aver$Date)

dygraph(CO7656_temp_xts) %>%
  dyAxis("y", label = "Daily temperature (°C)") 
#Average temperature by day for all water years:

CO7656_daily_wy_aver <- CO7656_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

CO7656_daily_wy_aver <- CO7656_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO7656_daily_wy_aver$aver_day_temp))

#str(CO7656_daily_wy_aver)
# try to show all years as means. 
CO7656_daily_wy_aver2 <- CO7656_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO7656_daily_wy_aver2$date_temp <- signif(CO7656_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO7656_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO7656_standard_dev <- CO7656_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO7656_standard_dev$residual)
## [1] 9.435791e-17

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO7656_standard_dev_all <- CO7656_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO7656_standard_dev_all <- CO7656_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7656_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.067411
1985 3.849404
1986 3.668068
1987 3.311225
1988 3.461834
1989 3.895538
1990 3.345777
1991 3.594486
1992 3.197390
1993 3.245250
1994 3.120369
1995 3.691098
1996 3.653801
1997 3.607704
1998 3.476696
1999 3.595767
2000 3.227173
2001 3.508776
2002 3.418585
2003 3.051783
2004 3.539072
2005 3.516482
2006 3.976319
2007 3.719573
2008 3.800601
2009 3.551654
2010 3.338182
2011 4.085668
2012 3.979993
2013 3.879710
2014 3.541149
2015 3.313898
2016 3.469157
2017 3.954068
2018 3.744968
2019 4.588709
2020 3.249373
2021 3.157533
2022 3.379058 
ggplot(CO7656_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7656_sd_mk <- mk.test(CO7656_standard_dev_all$sd_2)
print(CO7656_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO7656_standard_dev_all$sd_2
## z = 0.16936, n = 39, p-value = 0.8655
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 1.500000e+01 6.833667e+03 2.024291e-02
CO7656_sd_sens <- sens.slope(CO7656_standard_dev_all$sd_2)
print(CO7656_sd_sens)
## 
##  Sen's slope
## 
## data:  CO7656_standard_dev_all$sd_2
## z = 0.16936, n = 39, p-value = 0.8655
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.009483617  0.011403168
## sample estimates:
## Sen's slope 
## 0.001040035

Summer temperature standard deviation

CO7656_standard_dev_all_summer <- CO7656_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO7656_standard_dev_all_summer <- CO7656_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7656_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 1.935713
1985 2.119278
1986 2.013243
1987 2.073056
1988 2.498496
1989 1.833124
1990 2.237314
1991 1.703379
1992 2.008851
1993 1.979299
1994 1.776960
1995 2.618254
1996 1.742927
1997 1.796887
1998 2.076687
1999 1.851947
2000 1.491852
2001 2.012607
2002 1.815748
2003 1.581537
2004 1.951013
2005 1.782066
2006 3.519129
2007 1.753941
2008 1.807781
2009 1.796359
2010 1.912976
2011 1.552432
2012 3.558152
2013 1.794200
2014 1.579088
2015 1.982709
2016 2.006172
2017 1.680159
2018 3.389229
2019 4.043584
2020 2.006444
2021 1.973750
2022 1.710642 
ggplot(CO7656_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7656 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO7656_sd_mk_summer <- mk.test(CO7656_standard_dev_all_summer$sd_2)
print(CO7656_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO7656_standard_dev_all_summer$sd_2
## z = -0.84678, n = 39, p-value = 0.3971
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -71.00000000 6833.66666667   -0.09581646
CO7656_sd_sens_summer <- sens.slope(CO7656_standard_dev_all_summer$sd_2)
print(CO7656_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO7656_standard_dev_all_summer$sd_2
## z = -0.84678, n = 39, p-value = 0.3971
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.010744645  0.005620485
## sample estimates:
##  Sen's slope 
## -0.002716264

Winter temperature standard deviation

CO7656_standard_dev_all_winter <- CO7656_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO7656_standard_dev_all_winter <- CO7656_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO7656_standard_dev_all_winter <- CO7656_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7656_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.653469
1985 4.782512
1986 4.503770
1987 4.018840
1988 4.344342
1989 4.978024
1990 4.161851
1991 4.550557
1992 3.490449
1993 4.157226
1994 3.901893
1995 4.573380
1996 4.460319
1997 4.282828
1998 4.128093
1999 3.970478
2000 4.146942
2001 4.432268
2002 4.409185
2003 3.545802
2004 4.342837
2005 4.382139
2006 3.827248
2007 4.873215
2008 4.951051
2009 4.523490
2010 3.731486
2011 5.350289
2012 3.929436
2013 5.077884
2014 4.252081
2015 4.052701
2016 4.509374
2017 5.141089
2018 4.620516
2019 4.329715
2020 3.719626
2021 3.979182
2022 4.354394 
ggplot(CO7656_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7656 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO7656_sd_mk_winter <- mk.test(CO7656_standard_dev_all_winter$sd_2)
print(CO7656_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO7656_standard_dev_all_winter$sd_2
## z = -0.26613, n = 39, p-value = 0.7901
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -23.00000000 6833.66666667   -0.03103914
CO7656_sd_sens_winter <- sens.slope(CO7656_standard_dev_all_winter$sd_2)
print(CO7656_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO7656_standard_dev_all_winter$sd_2
## z = -0.26613, n = 39, p-value = 0.7901
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01337962  0.01320821
## sample estimates:
##  Sen's slope 
## -0.001632081

Spring and Fall temperature standard deviation

CO7656_standard_dev_all_spring <- CO7656_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO7656_standard_dev_all_spring <- CO7656_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7656_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.034073
1985 3.064609
1986 3.513850
1987 3.417461
1988 3.372144
1989 3.886313
1990 2.799211
1991 3.628112
1992 2.833636
1993 2.691351
1994 2.916113
1995 2.933359
1996 3.583731
1997 3.619193
1998 2.746210
1999 3.664649
2000 3.162393
2001 3.766521
2002 3.045279
2003 3.598301
2004 3.019071
2005 3.000378
2006 2.595234
2007 3.431935
2008 3.090242
2009 3.360528
2010 3.385341
2011 3.875988
2012 3.484504
2013 3.103193
2014 3.795843
2015 2.895605
2016 2.952351
2017 3.740011
2018 2.641638
2019 3.519195
2020 2.707509
2021 2.715416
2022 3.077034 
ggplot(CO7656_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7656 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO7656_sd_mk_spring <- mk.test(CO7656_standard_dev_all_spring$sd_2)
print(CO7656_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO7656_standard_dev_all_spring$sd_2
## z = -0.96775, n = 39, p-value = 0.3332
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  -81.0000000 6833.6666667   -0.1093117
CO7656_sd_sens_spring <- sens.slope(CO7656_standard_dev_all_spring$sd_2)
print(CO7656_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO7656_standard_dev_all_spring$sd_2
## z = -0.96775, n = 39, p-value = 0.3332
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.020572352  0.007441909
## sample estimates:
##  Sen's slope 
## -0.006201687

Fall temperature standard deviation

CO7656_standard_dev_all_fall <- CO7656_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO7656_standard_dev_all_fall <- CO7656_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO7656_standard_dev_all_fall <- CO7656_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7656_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.248861
1985 3.303691
1986 2.848955
1987 2.532503
1988 2.138187
1989 2.244281
1990 2.791397
1991 2.648730
1992 3.063703
1993 2.505885
1994 2.297185
1995 2.534642
1996 3.181316
1997 3.746296
1998 3.729575
1999 2.393939
2000 2.317650
2001 2.154905
2002 2.587524
2003 2.717277
2004 3.118313
2005 2.430619
2006 3.037807
2007 2.736879
2008 3.015199
2009 2.621915
2010 3.553546
2011 2.546956
2012 2.501661
2013 2.736491
2014 3.476960
2015 2.103975
2016 2.215745
2017 2.933247
2018 2.103867
2019 2.589082
2020 3.626047
2021 2.001434
2022 2.811885 
ggplot(CO7656_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7656 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO7656_sd_mk_fall <- mk.test(CO7656_standard_dev_all_fall$sd_2)
print(CO7656_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO7656_standard_dev_all_fall$sd_2
## z = -0.60484, n = 39, p-value = 0.5453
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -51.00000000 6833.66666667   -0.06882591
CO7656_sd_sens_fall <- sens.slope(CO7656_standard_dev_all_fall$sd_2)
print(CO7656_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO7656_standard_dev_all_fall$sd_2
## z = -0.60484, n = 39, p-value = 0.5453
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01757855  0.01131879
## sample estimates:
##  Sen's slope 
## -0.005088706

VALLECITO DAM CO8582

Detrending Data

#average water year temperature

CO8582_yearly_wy_aver <- COOP_stations_clean %>% 
  filter(station == "CO8582") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))

CO8582_temp_xts <- xts(CO8582_yearly_wy_aver$avg_T_c, order.by = CO8582_yearly_wy_aver$Date)

dygraph(CO8582_temp_xts) %>%
  dyAxis("y", label = "Daily temperature (°C)") 
#Average temperature by day for all water years:

CO8582_daily_wy_aver <- CO8582_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

CO8582_daily_wy_aver <- CO8582_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO8582_daily_wy_aver$aver_day_temp))

#str(CO8582_daily_wy_aver)
# try to show all years as means. 
CO8582_daily_wy_aver2 <- CO8582_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO8582_daily_wy_aver2$date_temp <- signif(CO8582_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO8582_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO8582_standard_dev <- CO8582_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO8582_standard_dev$residual)
## [1] 1.303179e-16

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO8582_standard_dev_all <- CO8582_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO8582_standard_dev_all <- CO8582_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8582_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.195055
1985 3.559954
1986 3.787658
1987 3.402581
1988 3.322995
1989 3.497396
1990 3.344103
1991 3.345133
1992 3.359204
1993 3.065636
1994 3.177175
1995 3.688527
1996 3.344606
1997 3.436496
1998 3.284087
1999 3.469742
2000 3.226386
2001 3.382710
2002 3.032005
2003 3.098054
2004 3.681084
2005 3.098605
2006 3.399274
2007 3.424611
2008 3.424672
2009 3.387154
2010 3.351483
2011 3.524069
2012 3.056001
2013 3.786870
2014 3.235817
2015 3.391219
2016 3.071576
2017 3.280991
2018 3.130957
2019 3.422241
2020 3.080326
2021 3.280592
2022 3.135074 
ggplot(CO8582_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8582_sd_mk <- mk.test(CO8582_standard_dev_all$sd_2)
print(CO8582_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO8582_standard_dev_all$sd_2
## z = -1.4516, n = 39, p-value = 0.1466
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -121.0000000 6833.6666667   -0.1632928
CO8582_sd_sens <- sens.slope(CO8582_standard_dev_all$sd_2)
print(CO8582_sd_sens)
## 
##  Sen's slope
## 
## data:  CO8582_standard_dev_all$sd_2
## z = -1.4516, n = 39, p-value = 0.1466
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.010288118  0.001644177
## sample estimates:
##  Sen's slope 
## -0.004135843

Summer temperature standard deviation

CO8582_standard_dev_all_summer <- CO8582_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO8582_standard_dev_all_summer <- CO8582_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8582_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 1.579077
1985 2.450251
1986 2.145256
1987 2.147308
1988 2.116747
1989 2.032895
1990 2.620783
1991 1.698953
1992 2.082496
1993 2.178575
1994 1.674607
1995 2.838381
1996 1.807302
1997 1.769535
1998 2.390738
1999 1.830583
2000 1.636308
2001 2.177360
2002 1.911096
2003 2.101649
2004 2.258015
2005 2.007679
2006 2.176184
2007 1.916562
2008 1.842907
2009 2.229087
2010 2.174356
2011 1.572073
2012 1.545226
2013 1.834963
2014 1.793817
2015 2.065234
2016 2.287603
2017 1.886518
2018 1.606058
2019 2.206086
2020 2.012836
2021 2.068044
2022 2.257185 
ggplot(CO8582_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8582 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO8582_sd_mk_summer <- mk.test(CO8582_standard_dev_all_summer$sd_2)
print(CO8582_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO8582_standard_dev_all_summer$sd_2
## z = -0.48387, n = 39, p-value = 0.6285
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -41.00000000 6833.66666667   -0.05533063
CO8582_sd_sens_summer <- sens.slope(CO8582_standard_dev_all_summer$sd_2)
print(CO8582_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO8582_standard_dev_all_summer$sd_2
## z = -0.48387, n = 39, p-value = 0.6285
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.011667762  0.005849158
## sample estimates:
##  Sen's slope 
## -0.001907898

Winter temperature standard deviation

CO8582_standard_dev_all_winter <- CO8582_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO8582_standard_dev_all_winter <- CO8582_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO8582_standard_dev_all_winter <- CO8582_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8582_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.847057
1985 4.049547
1986 4.494387
1987 4.256319
1988 3.821482
1989 4.442285
1990 4.094723
1991 4.358364
1992 3.749949
1993 3.747999
1994 3.926086
1995 4.321029
1996 3.987300
1997 4.038893
1998 3.756625
1999 3.572155
2000 4.126135
2001 3.981876
2002 3.492047
2003 3.529586
2004 4.485568
2005 3.849973
2006 4.135002
2007 4.304422
2008 4.321792
2009 4.048524
2010 3.631638
2011 4.436681
2012 3.769156
2013 4.877975
2014 3.693727
2015 3.889881
2016 3.696468
2017 3.853385
2018 3.984933
2019 3.876716
2020 3.222906
2021 3.914995
2022 3.775318 
ggplot(CO8582_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8582 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO8582_sd_mk_winter <- mk.test(CO8582_standard_dev_all_winter$sd_2)
print(CO8582_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO8582_standard_dev_all_winter$sd_2
## z = -1.4274, n = 39, p-value = 0.1535
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -119.0000000 6833.6666667   -0.1605938
CO8582_sd_sens_winter <- sens.slope(CO8582_standard_dev_all_winter$sd_2)
print(CO8582_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO8582_standard_dev_all_winter$sd_2
## z = -1.4274, n = 39, p-value = 0.1535
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.017213348  0.001836151
## sample estimates:
##  Sen's slope 
## -0.006983808

Spring and Fall temperature standard deviation

CO8582_standard_dev_all_spring <- CO8582_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO8582_standard_dev_all_spring <- CO8582_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8582_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.705329
1985 3.138231
1986 3.395151
1987 2.932670
1988 3.849246
1989 3.133287
1990 2.437669
1991 3.118406
1992 2.821949
1993 2.499646
1994 3.351265
1995 2.995433
1996 3.488868
1997 3.460010
1998 3.073371
1999 3.602936
2000 3.086724
2001 3.526095
2002 2.738884
2003 3.504942
2004 2.833372
2005 3.128759
2006 2.735151
2007 3.397614
2008 2.851106
2009 3.384965
2010 3.238152
2011 3.660475
2012 3.151708
2013 3.343300
2014 3.626165
2015 3.091006
2016 2.667774
2017 3.159069
2018 2.906205
2019 3.556864
2020 2.598872
2021 2.951218
2022 2.544456 
ggplot(CO8582_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8582 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO8582_sd_mk_spring <- mk.test(CO8582_standard_dev_all_spring$sd_2)
print(CO8582_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO8582_standard_dev_all_spring$sd_2
## z = -0.89517, n = 39, p-value = 0.3707
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  -75.0000000 6833.6666667   -0.1012146
CO8582_sd_sens_spring <- sens.slope(CO8582_standard_dev_all_spring$sd_2)
print(CO8582_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO8582_standard_dev_all_spring$sd_2
## z = -0.89517, n = 39, p-value = 0.3707
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.017914984  0.006302449
## sample estimates:
##  Sen's slope 
## -0.006452656

Fall temperature standard deviation

CO8582_standard_dev_all_fall <- CO8582_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO8582_standard_dev_all_fall <- CO8582_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO8582_standard_dev_all_fall <- CO8582_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8582_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.369152
1985 3.222177
1986 3.221453
1987 2.869632
1988 2.806946
1989 2.259131
1990 2.602142
1991 2.405738
1992 3.429397
1993 2.755296
1994 2.525678
1995 3.201853
1996 2.917753
1997 3.705080
1998 3.148361
1999 1.986162
2000 2.582096
2001 2.425238
2002 2.613933
2003 2.446174
2004 3.713580
2005 2.214157
2006 3.085811
2007 2.619196
2008 2.724464
2009 2.882420
2010 3.475592
2011 2.410535
2012 2.333219
2013 2.556643
2014 2.837810
2015 2.276561
2016 2.505294
2017 2.749786
2018 2.461669
2019 2.873154
2020 3.855550
2021 2.638235
2022 3.121268 
ggplot(CO8582_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8582 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO8582_sd_mk_fall <- mk.test(CO8582_standard_dev_all_fall$sd_2)
print(CO8582_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO8582_standard_dev_all_fall$sd_2
## z = -0.096775, n = 39, p-value = 0.9229
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##   -9.00000000 6833.66666667   -0.01214575
CO8582_sd_sens_fall <- sens.slope(CO8582_standard_dev_all_fall$sd_2)
print(CO8582_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO8582_standard_dev_all_fall$sd_2
## z = -0.096775, n = 39, p-value = 0.9229
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01413440  0.01239632
## sample estimates:
##  Sen's slope 
## -0.001128856

DEL NORTE 2E CO2184

Detrending Data

#average water year temperature

CO2184_yearly_wy_aver <- COOP_stations_clean %>%
  filter(station == "CO2184") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO2184_daily_wy_aver <- CO2184_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(avg_T_c))

#average mean temperature by day for the period of record:

CO2184_daily_wy_aver <- CO2184_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO2184_daily_wy_aver$aver_day_temp))

#str(CO2184_daily_wy_aver)
# try to show all years as means. 
CO2184_daily_wy_aver2 <- CO2184_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO2184_daily_wy_aver2$date_temp <- signif(CO2184_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO2184_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO2184_standard_dev <- CO2184_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO2184_standard_dev$residual)
## [1] -2.097473e-16

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO2184_standard_dev_all <- CO2184_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO2184_standard_dev_all <- CO2184_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO2184_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.798354
1985 3.238769
1986 3.730457
1987 3.067875
1988 3.258441
1989 3.537167
1990 3.106260
1991 3.338810
1992 3.473238
1993 3.065377
1994 3.249189
1995 3.714892
1996 3.667398
1997 3.612252
1998 2.940217
1999 3.685857
2000 3.341569
2001 3.096662
2002 3.105707
2003 3.308729
2004 3.736161
2005 3.426736
2006 3.472574
2007 3.707045
2008 3.856440
2009 3.158845
2010 3.311903
2011 3.701188
2012 3.346574
2013 4.305294
2014 3.358905
2015 3.343961
2016 3.151279
2017 3.403460
2018 3.099286
2019 3.223172
2020 3.436846
2021 3.341378
2022 3.480513 
#CALLING THIS something different
CO2184_all_V2 <- ggplot(CO2184_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')


print(CO2184_all_V2)

CO2184_sd_mk <- mk.test(CO2184_standard_dev_all$sd_2)
print(CO2184_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO2184_standard_dev_all$sd_2
## z = 0.24194, n = 39, p-value = 0.8088
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 2.100000e+01 6.833667e+03 2.834008e-02
CO2184_sd_sens <- sens.slope(CO2184_standard_dev_all$sd_2)
print(CO2184_sd_sens)
## 
##  Sen's slope
## 
## data:  CO2184_standard_dev_all$sd_2
## z = 0.24194, n = 39, p-value = 0.8088
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.007626337  0.007591167
## sample estimates:
##  Sen's slope 
## 0.0004961876

Summer temperature standard deviation

CO2184_standard_dev_all_summer <- CO2184_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO2184_standard_dev_all_summer <- CO2184_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO2184_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 1.701728
1985 2.069629
1986 1.675529
1987 1.686232
1988 1.652258
1989 1.564582
1990 2.116643
1991 1.522834
1992 1.885193
1993 1.835922
1994 1.569129
1995 2.292973
1996 1.541520
1997 1.739051
1998 1.993972
1999 1.749043
2000 1.630222
2001 1.656631
2002 1.692479
2003 2.130000
2004 2.195237
2005 1.852805
2006 1.959780
2007 1.732923
2008 1.877601
2009 1.785384
2010 1.874888
2011 1.476163
2012 1.492401
2013 2.184176
2014 1.864381
2015 1.769342
2016 1.978573
2017 1.924525
2018 1.758846
2019 1.835296
2020 1.858668
2021 2.310526
2022 2.024805 
ggplot(CO2184_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO2184 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO2184_sd_mk_summer <- mk.test(CO2184_standard_dev_all_summer$sd_2)
print(CO2184_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO2184_standard_dev_all_summer$sd_2
## z = 1.7419, n = 39, p-value = 0.08152
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  145.0000000 6833.6666667    0.1956815
CO2184_sd_sens_summer <- sens.slope(CO2184_standard_dev_all_summer$sd_2)
print(CO2184_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO2184_standard_dev_all_summer$sd_2
## z = 1.7419, n = 39, p-value = 0.08152
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.001211471  0.011159386
## sample estimates:
## Sen's slope 
## 0.005774722

Winter temperature standard deviation

CO2184_standard_dev_all_winter <- CO2184_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO2184_standard_dev_all_winter <- CO2184_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO2184_standard_dev_all_winter <- CO2184_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO2184_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.286967
1985 3.823302
1986 4.731115
1987 3.853111
1988 3.862636
1989 4.395363
1990 3.576006
1991 4.383930
1992 3.373851
1993 3.669747
1994 3.888883
1995 4.195896
1996 4.300951
1997 4.304595
1998 3.519306
1999 3.852622
2000 3.841239
2001 3.879165
2002 3.777837
2003 3.972689
2004 4.593027
2005 4.382614
2006 4.311472
2007 4.840707
2008 4.772901
2009 3.863896
2010 3.654145
2011 4.758157
2012 4.042888
2013 5.174314
2014 4.009521
2015 4.058873
2016 3.774087
2017 3.996472
2018 3.642315
2019 3.876825
2020 3.692999
2021 3.535686
2022 4.390758 
ggplot(CO2184_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO2184 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO2184_sd_mk_winter <- mk.test(CO2184_standard_dev_all_winter$sd_2)
print(CO2184_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO2184_standard_dev_all_winter$sd_2
## z = 0.072581, n = 39, p-value = 0.9421
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 7.000000e+00 6.833667e+03 9.446694e-03
CO2184_sd_sens_winter <- sens.slope(CO2184_standard_dev_all_winter$sd_2)
print(CO2184_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO2184_standard_dev_all_winter$sd_2
## z = 0.072581, n = 39, p-value = 0.9421
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01294191  0.01123660
## sample estimates:
##  Sen's slope 
## 0.0004791036

Spring and Fall temperature standard deviation

CO2184_standard_dev_all_spring <- CO2184_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO2184_standard_dev_all_spring <- CO2184_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO2184_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.216438
1985 2.678631
1986 2.997311
1987 3.046403
1988 3.248865
1989 3.510275
1990 2.710673
1991 2.998136
1992 2.724351
1993 2.427897
1994 3.301673
1995 3.031827
1996 3.720276
1997 3.522047
1998 2.506387
1999 3.481011
2000 3.752910
2001 2.881770
2002 2.666107
2003 3.562739
2004 3.208543
2005 3.134278
2006 2.940810
2007 3.523360
2008 3.268286
2009 2.967166
2010 3.502197
2011 3.660154
2012 3.341963
2013 3.794801
2014 3.641980
2015 2.809653
2016 3.299330
2017 3.491309
2018 2.449094
2019 3.479622
2020 2.618257
2021 3.189132
2022 3.273935 
ggplot(CO2184_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO2184 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO2184_sd_mk_spring <- mk.test(CO2184_standard_dev_all_spring$sd_2)
print(CO2184_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO2184_standard_dev_all_spring$sd_2
## z = 0.62904, n = 39, p-value = 0.5293
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 5.300000e+01 6.833667e+03 7.152497e-02
CO2184_sd_sens_spring <- sens.slope(CO2184_standard_dev_all_spring$sd_2)
print(CO2184_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO2184_standard_dev_all_spring$sd_2
## z = 0.62904, n = 39, p-value = 0.5293
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.008274262  0.017352568
## sample estimates:
## Sen's slope 
## 0.005268223

Fall temperature standard deviation

CO2184_standard_dev_all_fall <- CO2184_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO2184_standard_dev_all_fall <- CO2184_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO2184_standard_dev_all_fall <- CO2184_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO2184_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 1.994998
1985 2.952172
1986 2.314373
1987 2.442785
1988 2.404126
1989 2.590376
1990 2.716540
1991 2.637648
1992 3.122279
1993 2.819353
1994 3.249444
1995 2.411606
1996 2.933987
1997 3.539284
1998 2.932876
1999 2.618765
2000 2.854073
2001 2.567835
2002 2.682918
2003 2.198283
2004 3.530130
2005 2.468374
2006 3.002252
2007 2.562223
2008 2.777312
2009 2.836104
2010 3.581238
2011 2.785499
2012 2.563709
2013 2.991718
2014 2.805952
2015 2.354842
2016 2.172908
2017 2.820517
2018 2.728988
2019 2.303466
2020 4.021269
2021 3.721992
2022 2.877937 
ggplot(CO2184_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO2184 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO2184_sd_mk_fall <- mk.test(CO2184_standard_dev_all_fall$sd_2)
print(CO2184_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO2184_standard_dev_all_fall$sd_2
## z = 1.2339, n = 39, p-value = 0.2172
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  103.0000000 6833.6666667    0.1390013
CO2184_sd_sens_fall <- sens.slope(CO2184_standard_dev_all_fall$sd_2)
print(CO2184_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO2184_standard_dev_all_fall$sd_2
## z = 1.2339, n = 39, p-value = 0.2172
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.00389730  0.02075506
## sample estimates:
## Sen's slope 
## 0.008219316